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Abstract 

Recent observations in X-ray and Gamma-ray suggest that the emission region of 
the pulsar magnetosphere can be multifold. In particular, the open-close boundary of 
the magnetic field, so-called the Y-point, can be the place where magnetic field energy 
converts into the plasma heat and/or flow energy. Here, we present a new Particle- 
in-Cell code in axisymmetric geometry, which can be applied to the Y-point of the 
pulsar magnetosphere. The electromagnetic solver is used in the two-dimensional grid 
points with the cylindrical coordinate (-R, z), while the particle solver operates in the 
three-dimensional Cartesian coordinate (x, z) by use of Buneman-Boris method. 
The particle motion can be relativistic. The inner boundary conditions are set up to 
generate rotation of the magnetosphere in use of the force-free semi-analytic solution 
given by Uzdensky (2003, ApJ, 598, 446). 

The code has been verified by dispersion relations of all the wave modes in an 
electron-positron plasma. 

The initial test run is also presented to demonstrate the Y-shaped structure at 
the top of the dead zone on the light cylinder. We suggest that the structure is 
variable with quasi-periodicity with magnetic reconnection and that plasma will be 
accelerated and/or heated. In the time-averaged point of view, break up of the ideal- 
MHD condition takes place in the vicinity of the Y-point. 
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1. Introduction 

Recently many authors have intensively studied the magnetosphere of the rotation pow- 
ered pulsars by means of computer simulations in force- free approximation (Contopoulos et 
al. 1999; Mckinney 2006; Timokhin 2006; Spitkovsky 2006), magneto-hydrodynamics (MHD) 
(Komissarov 2006; Bucciantini et al. 2006), two component MHD (Kojima & Oogi 2009) and 
particle simulation (Wada & Shibata 2007). In such studies, boundary layers, e.g. the equa- 
torial current sheet and the surface of the star, are sometimes troublesome. Dissipation is 
suggested in the boundary layers (Komissarov 2006), but it can be numerical. Acceleration 
and heating likely take place in the magnetic neutral sheet around the equatorial plane. The 
inner edge of the neutral sheet, which is called the "Y-point", is the open-close boundary of 
the magnetic field structure (figure 1). Location of the Y-point is not certain and can be inside 
of the hght cyhnder (Uzdensky 2003). 

In this paper, we study the Y-point in microscopic point of view to understand the 
dissipative process in the boundary layer. We expect that our simulation provides a hint 
for the appropriate boundary conditions for the global simulation. As has been suggested by 
Lyubarskii (1996), another possibility is that the Y-point might be a gamma-ray source. Pulsed 
gamma-ray radiation is previously attributed to the outer gaps and/or the polar caps. One 
may say that the radiation near and beyond the light cylinder becomes DC components, but 
this may not be the case. Recently, Kirk et al. (2002) proposed that the magnetic neutral sheet 
in the wind far from the light cylinder can be the source of the gamma-ray pulses. If dissipation 
of the magnetic field around the Y-point takes place, radiation from the vicinity of the Y-point 
would also contribute to the pulsed high energy emission. 

Uzdensky (2003) calculated the force-free solution around the Y-point. However, he 
found that the force- free condition is broken around the neutral sheet. The global MHD simu- 
lation (Komissarov 2006) also suggested breakdown of force-free condition around the neutral 
sheet and the magnetic energy conversion to the thermal and kinetic energy. The plasma inertia 
and magnetic dissipation seem to play an important role around the Y-point. 

We intend to study the Y-point via newly developed axisymmetric particle-in-cell (PIC) 
method in the cylindrical coordinates. This paper is aimed to establish the numerical code, the 
results of code check, and the initialization for the Y-point simulation. Some of the initial result 
shall be given though detailed analysis shall be given in a subsequent paper. In our scheme, 
the system is reflection symmetry with respect to the equator, which is the lower boundary. 
The free boundary condition is imposed on the outer boundary, which should be continued to 
the wind zone. The left and upper boundaries (inner boundaries) are assumed to satisfy the 
force-free solution given by Uzdensky (2003). 
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Fig. 1. The meridional plane of the pulsar magnetosphere. 
2. Particle-in-Cell Code 

Both particle code and Vlasov code describe the plasma kinetic effect. However, Vlasov 
simulation requires huge computer memory, so that we use the particle code in our simulation. 
The simplest system to treat the Y-point would be the axisymmetric one, and therefore we 
developed axisymmetric PIC code using cylindrical coordinates. 

2.1. Particle motion solver 

We use Buneman-Boris method in Cartesian coordinates. Then, particles are distributed 
in the three-dimensional Cartesian coordinates {x,y,z) and solved by Buneman-Boris method. 
Therefore, if a suitable inner boundary condition of E and B are given, particles rotate around 
the z-a.xis in the 3D space. Afterwords, we transform the particle coordinates to the cylindrical 
ones {R, (f, z) and solve the field quantities in the R-z plane. 
The relativistic equation of motion is 
du q 
dt m 

where u is the space part of the four-velocity, q is the charge, m is the rest mass, 7 is the 
Lorentz factor, and E and B are the electric and the magnetic fields, respectively. In Cartesian 
coordinate, E and B are given by 

E = {Ejicosip — E^smip,ERsiinp + E^cos(p,Ez), (2) 
B = {BjiCos(p — B^pSm(p,Bjismip + B^pCos(p,Bz). (3) 

Equation (1) is rewritten in a centered-difference form as follows 



E-\ u X B 

C7 



(1) 
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= ^ Ie" + ^ X , (4) 

At m V 2c7 ; ' ^ ^ 

where n and At are the number and increment of the time steps. 

Buneman-Boris method divides equation (4) into three steps. 
1st step : an acceleration by E for At/2 

m 2 ^ ^ 

2nd step : a rotation by B for At 

u^ = u~ + u~ X T, (6) 

u^ = u- + -^^u' X T, (7) 



where T = qB''At/{2mc-f-) and 7- = y'l + (n-/c)2. 
3rd step : an acceleration by E for At/2 

u^+y^ = u+ + ^E^—. (8) 
m 2 ^ ^ 

Finally, 1;"+^/^ is calculated as follows, 

^n+1/2 _ t,"+l/Yy^l + (^n+l/2/c)2. (9) 

The position of the particles is updated by 

^n+l^^„+l/2^^n+l/2At^ (^q) 

Electromagnetic field solver 

Maxwell equations are solved in the cylindrical coordinates with axisymmetry by the 
Leap-Frog method. The magnetic field is updated by the Faraday's law, 
f) R 

-^-cVxE (11) 

and the electric field is updated by the Ampere's law, 
f) V 

— = -AttJ + cVxB, (12) 
ot 

where J is the current density. 

2.2.1. Field definition 

In the R-z plane, fields are defined as shown in figure 2, with indices j and k for R 
and z coordinates, respectively. For example, the (/^-components of E, B and J are given at 
{Rj+i/2, Zk+1/2), but the /^-components are given at {Rj+i/2,Zk). 

2.2.2. Update of magnetic field 

We update the magnetic field by At/2 using the Faraday's law. 
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Fig. 2. Location of the grid quantities. 

dt ^[ dz dRj' 



dt RdR'" 
which are written in centered-difference form, 

^R,j+l/2,k — ^R,j+l/2,k 

^^,j+l/2,k+l/2 ~ ^^,j+l/2,k-l/2 At 



+c- 



Az 2 



^ip,j+l/2,k+l/2 - ^^,j+l/2,k+l/2 

Crpn rpn rpn rpn '■ 

-'^i?j+l/2,fc+l -^R,j+l/2,k ^z,j+l,k+l/2 ^z,j,k+l/2 

C -^i+l/2-5'.^,j+l/2,fc+l/2 ~ -^j-l/2-^y,j-l/2,fc+l/2 At 

'Rj 'Kr 2"' 

^.^.5'. Update of electric field 

We update the electric field by At using the Ampere's law, 

^=-'"''^-^^' 



dt ^ \ dz OR J' 

^ = -4TrJ, + -^(RBj, 
dt RdR^ 

which are rewritten in centered-difference form, 

pn+l pin 

^R,j+l/2,k — ^R,j+l/2,k 



T^n+1/2 _ nn+1/2 

. jn+1/2 ^ -Py,j+l/2,fc+l/2 -Py.i+l/2,fc-l/2 , . 



pn+1 pn 

-^(^,j+l/2,fc+l/2 ~ -^(/3j+l/2,A:+l/2 



_4 7n+l/2 , 
^^'^V',i+l/2,fc+l/2 ^^ c 



^z,j,k+l/2 — ^z,j,k+l/2 



B 



n+l/2 



B 



n+l/2 



R,j + l/2,fc+l -"flj+l/2,fc 



B 



n+l/2 



B 



n+l/2 



z,j+ l,k+l/2 ^z,j,k+l/2 

AR 



R R"+V2 p Rn+1/2 

_4^jn+l/2 ^ _^ C -^i+l/2^yj+i/2,fc+l/2 -^i-l/2-Py,j-l/2,fc+l/2 | 



AR 



(22) 

At, 
(23) 

(24) 



2.3. Calculation of charge density 

To obtain the charge density at {Rj,Zk), charge of a particle is distributed to the neigh- 
boring four grids with weightings in proportion to the volumes opposite to the particle, i.e. the 
shaded region in figure 3 for the left-bottom grid point. We calculate thus the charge density 




particle 



R 



j+i 



Fig. 3. Schematic diagram for interpolation between a particle and grids, 
with shape factors. 



Sj {Ri) 



{Rj+i - Ri) /AR for Rj < Ri < Rj+i 
{Ri - Rj^i) /AR for Rj^i < Ri < Rj 
otherwise 



{{zk+i - Zi)/Az for Zk <Zi< Zk+i 
{zi - Zk-i) /Az for Zk-i < Zi< Zk 
otherwise 
in our code. The charge distributed to {Rj,Zk) is 

Qj,k = ^QiSj (Ri) Sk (Zi) , 

i 
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(25) 



(26) 



(27) 



which yields the charge density, 

Pj,k = Qj,k/Vj^k, (28) 

where V,- ^ = 27rRjARAz. 

On the other hand, the electric field at the position of the i-th particle {Ri,Zi) is given 

by 



E{Ri,Zi) — SjSk 



Sj+iSk 



( {Rj+i/Ri)ERj+i,k+i \ 
Sj+iSk+i E^j^i^k+1 



J 



Ez,j+l,k 

The magnetic field is also given in the same form. 
2.4- Calculation of current density 

Current density is also calculated with the shape factors: 

'^hk = J^^iSj (Ri) Sk (Zi) Vi. 



\ 



E. 



(29) 



z,j+l,k+l 



J 



(30) 



We interpolate the current densities at the full integer grids to those at the appropriate 
grid points shown in figure 2 by 

Rj+lJ R,j+l,k + RjJR,j,k 



jR,j+l/2,k 
Jip,j+l/2,k+l/2 



J. 



z,j,k+l/2 



2-Rj+i/2 

Jip,j+l,k+l ~l~ Jip,j,k+1 ~l~ Jip,j-\-l,k ~l~ Jtp,j,k 

4 

Jz,j,k+1 ~l~ Jz,j,k 



(31) 

(32) 
(33) 



Since the above current density does not guarantee the charge conservation law, the 
electric fields contain an error. Therefore, we use the following procedure to correct the electric 
field. Gauss' law yields 

V-{E' + E,)=47rp, (34) 

where E' is the obtained field by (22), (23), (24) and Ec is the correction. It follows from (34) 
that Ec is obtained by 

V-E, = 47rpc, (35) 
where Pc = p — V ■ E' /An. Then, we solve 

-V^c = 47rpe (36) 
and we get E^ = — V0c- 
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2. 5. A calculation of the PIC method for one time step 

A calculation cycle of the PIC method is shown in figure 4. We interpolate the electro- 
magnetic fields at the particle positions, -Ej, Bi, from the values of the given grids points, Ej^^, 
Bj^k- We integrate the equations of motion and move the particle positions. We next calculate 
the charge densities and current densities using the particle data. Electromagnetic fields are 
updated by integrating Maxwell equations. 



Integration of equations of motion 
and moving particles 




Calculation of p, J 



Integration of Maxwell equations 



Fig. 4. A calculation cycle of the PIC method. 



3. Result 

3.1. Code check 

In order to check the code, we obtain the dispersion relations of waves in electron- 
positron plasmas for Light mode, X-mode, Fast mode, 0-mode, Electromagnetic mode, and 
Alfven mode. The boundary conditions are perfect conductor for the outer radius and periodic 
for ^-direction. The plasma is uniform and has no bulk velocity. Particle distribution for the 
thermal velocities is non-relativistic Maxwellian. The simulation settings are given in table 1. 
3.1.1. Light mode 

The simulation for Case A in Table 1 demonstrates the light-mode. 

The linear dispersion relation of the light mode is 

uj^ = ujl + ^k\ (37) 



where 



where 



m7th 



(38) 



7th = l + T^ — ^, 39 

i — 1 mc^ 

r is a ratio of specific heat, is the Boltzmann constant, and T is the temperature. The 
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Table 1. The simulation settings for the code check. Bo is background magnetic field, iith is thermal velocity, ujp is plasma 
frequency, and SI is gyro-frequency. 





Case A 


Case B 


Case C 


Bo 


(0, 0, 0) 


(0, B^, 0) 


(0, 0, B^) 


Vth 


0.1c 


0.1c 


0.1c 


AR,Az 


1 


1 

VthOJp 


1 

VthOJp 


At/AR 


0.2 


0.2 


0.2 


Grid number {R x z) 


128 X 512 


128 X 512 (X,Fast) 
256 X 512 (0) 


128 X 512 


Total time step 


4096 


4096 


4096 


Particle number per grid 


200 


200 


200 


Cl/ojp 





1 


1 



frequency of the light mode is cUp at A; = and asymptotically it approaches to cu = ck in high 
wave numbers. 

We compute the Fourier transform with respect to z and then integrate that in in- 
direction. We do not impose any kind of initial perturbation, i.e., there is only thermal noise 
initially. Figure 5 shows the result for B^, showing that dispersion relation (37) is well repro- 
duced. 

3.1.2. X-mode, Fast mode, and 0-mode 

We next analyze the perpendicular propagating waves. The simulation settings are given 
in Case B in Table 1. 

The linear dispersion relations of X-mode (positive-signed) and Fast mode (negative- 
signed) are 



2 1 

" =2 



c' + ci)k' + n'+uj'p±J{c^- cifk^ + 2(c2 - cI){ujI -n^)k^ + {n^ + ujjy ,(4o) 



where 



eB 

n = ^^, (41) 



mc7th 

and Cg is the speed of the sonic wave, 



/(r-l)2:H^c. (42) 

7th 



The frequency of the X-mode becomes the upper hybrid frequency cuuh = y^p + at 
k = 0, and asymptotically approaches ui = ck in high wave numbers. The frequency of the Fast 
mode is a; = at A; = 0, and asymptotically approaches to u = Cgk in high wave numbers. 

Figure 6 shows Fourier transform of B^^, showing that these waves are well fitted by the 
dispersion relations (40). 
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The linear dispersion relation of 0-mode is 



2 2 I 2/2 

u = uj^ + c k . 



(43) 



The frequency of the 0-mode is cjp at A; = 0, and asymptotically approaches to oo = ck 
for high wave numbers. 

Figure 7 shows Fourier transform of showing that this wave is well fitted by the 
dispersion relation (43). 

3.1.3. Electromagnetic mode and Alfven mode 

Finally, we analyze the parallel propagating waves, Case C, for which the simulation 
settings are given in Table 1. 

The linear dispersion relations of the electromagnetic mode (positive-signed) and Alfven 
mode (negative-signed) are 



The frequency of the electromagnetic mode is cjuh at k = 0, and asymptotically ap- 
proaches to uj = ck in high wave numbers. The frequency of the Alfven mode is = at /c = 0, 
and asymptotically approaches to the relativistic gyro-frequency in high wave numbers. 

Figure 8 shows Fourier transform of B^, showing that these waves are well fitted by the 
dispersion relations (44). 
3.1.4- Summary of the code check 

We perform the simulations for waves in electron-positron plasmas in order to check 
the newly developed axisymmetric PIC code. We find that dispersion relations of Light mode, 
X-mode, Fast mode, 0-mode, Electromagnetic mode, and Alfven mode are well reproduced. 
Thus we conclude that the code is successfully constructed. 

3.2. Application to the Y-point of the pulsar magnetosphere 

3.2.1. Outline of the Y-point simulation 

Regarding the local simulation in the global structure shown in figure 1, the boundary 
condition is of great importance for the simulation of the Y-point. Inner region would be 
well approximated by the force-free solution, so that the left and upper boundaries (inner 
boundaries) are assumed to be the force-free solution by Uzdensky (2003) as shown in figure 
9. On the other hand, it is known that Uzdensky 's solution is broken down beyond the light 
cylinder. It is likely that the particle inertia becomes important, and the frozen- in condition 
breaks down near the Y-point where the magnetic energy density decreases and plasma is 
more important. The outer (right-side) boundary may be crossed by super-fast wind. Then 
we impose free boundary for the right-side boundary. The lower boundary can be treated as 
reflection symmetry. The boundary conditions might be fair only if the simulation box is large 
enough. As for the initial condition we also use the Uzdensky's solution. The simulation shall 
be continued until the system becomes quasi steady state. 




(44) 
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Fig. 5. The dispersion relation diagrams of tlic light mode. Gray map indicates Fourier transformation 
of B^. The left figure is the full simulation scale and the right is the lower frequency part. 
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10 12 14 



1 
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p 



Fig. 6. Fourier transformation of B^p. The dispersion relation diagrams of the X-mode (upper frequency 
wave) and the Fast mode (lower frequency wave). The left figure is the full simulation scale and the right 
is the lower frequency part. For this case, wuh = and Cg ~ 6.6 x 10~^. 

3.2.2. The force-free solution around the Y-point 

As shown by Uzdensky (2003), if the Y-point is located at the hght cyhnder, the regular- 
ity condition at the light cylinder determines the poloidal magnetic flux , and the poloidal 
current, /(\&). The obtained magnetic field is shown in figure 10. Once the force-free solution 
\l/ and /(\1/) are obtained, the electromagnetic field, charge density, and current density are 
given by 
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Fig. 8. Fourier transformation of B^p for Case C. The dispersion relation of the Electromagnetic mode 
(upper frequency wave) and the Alfven mode (lower frequency wave) can be seen. The left figure is the 
full simulation scale and the right is the lower frequency part. 

Bp = X Vv^, (45) 

= |, (46) 

E = -'^V^, (47) 
c 

p=i-V-^;, (48) 
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(Field boundary condition) 



(Particle boundary condition) 



force free 



remove or inject so that p. ^ and X ^ 
keep force free 




5 / SR = 



equator 



symmetry 




remove 



equator 



Fig. 9. The boundary conditions for the Y-point simulation. 




Fig. 10. The force-free solution around the Y-point given by Uzdensky (2003). The dashed lines represent 
poloidal magnetic field line. The solid line represents the separatrix. 'J = 0, between the open and closed 
field-line regions. 



c 



(49) 
(50) 



where VL is the angular velocity of the star. 

The bulk velocities of the plasma are restricted by 

i;± = r2 X r (51) 

where -|- and - respectively indicate positron and electron, r is the position vector, Q, is the 
angular velocity vector, and k is a scalar function to be determined. 
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Since the force-free solution does not concern about the plasma density and bulk velocity, 
we introduce two model parameters, 

n = n+ + n_ (52) 

and 

1= I ^ (53) 

where v = |(7+n+v+ + 7_n_'y_)/ (7+n+ + 7_n_)|. Giving the two model parameters, we obtain 
n± and v±, which satisfy the force-free solution for p and J. The Nj^k = n±yj,k particles, whose 
velocity is v±, are distributed uniformly in the (j,A;)-cell by a random number generator. 
3.2.3. Inner boundary condition 

The inner boundary conditions are imposed for the innermost and uppermost cells (figure 
11). For the inner boundary cells, we assume the density, n, so that the number of particle in a 
cell, iVj = nVj^ki is constant. With the force-free charge density, we have the particle numbers 
for each species to be distributed in each cell, 

N^ = \{n,,u±^V,,^, (54) 

where we typically take AT,- fc=constant=100. The field-aligned velocity of the particles are 
restricted by the force- free current density, i.e., 

g(n+«:+-n„«:„) = |-^, (55) 
where the poloidal velocities for each species are given by 

= K±Bp. (56) 

The mean value of k is given by 

_ ^ 7+ra+/t+ + 7_n_ft:_ 
7+n+ + 7_n_ 



where 7-1- = l — {v±/c)'^ ' . Given 7, (55) and (57) are the simultaneous equations with 
respect to k±. In the present simulation, we set typically 7 = 50. Thus the velocity in the 
boundary cells are given by (51). The current density J = q{n^v^ — n^vJ) automatically 
satisfies both (49) and (50). 

We remove or inject particles in the innermost and uppermost cells so that p± and J± 
keep these force-free formulae. 
3.2.4. Initialization 

As was mentioned before, the initial state of the plasma is set to be the Uzdensky's 
force-free solution. The electromagnetic field is set by (45), (46) and (47). For the particles, 
we defined the closed field region, i.e. the dead zone, as \l/ > and the open field region, i.e. 
the wind zone, as \l/ < 0. In the dead zone, we take A'j-^fc=100 and k±=0. It is known that there 
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-1/2 




Fig. 11. The innermost cells and the uppermost cells. We imposed force-free condition for those cells. In 
the wedge shape region 6 <9c, Uzdensky solution breaks down. 

is a critical line shown in figure 11 with the critical angle 9c ^ 62° from the equatorial plane 
and that in the wedge like region 9 < 9c, we have l^^l > |.B|. In the force-free wind zone, we 
use the continuity condition, V ■ {nj-v±) = 0, instead of setting A'^-^^, and take 7 = 50. In the 
non-force-free wind zone, we assume Vp± \\ Bp instead of using (51) , and take Nj k=^00 and 7 
= 50. In any case, the memory of the initial condition will be lost in a few light-transit times, 
so that the final state will not be affected by the initial condition. 
3.2.5. Normalization 

We normalize the distance, velocity, and magnetic field respectively by the light radius, 
i?LC = c/fi, the speed of light, and (i?*/2)(i?*/i?Lc)^! where is the field strength at the poles 
and i?^, is the radius of the star, respectively. In order to increase computational efficiency, we 
renormalize distance and time by the grid spacing. Ax, and At/2, respectively. 

Normalization constants and renormalization constants are shown in table 2. 

3.3. Initial result of the Y-point simulation 

We have performed a test run for the Y-point. After a time of 1/60 rotation, the plasma 
attains a steady state, in some sense, with quasi-periodic variability with time scale of 1/300 
rotation. Figure 12 shows the poloidal magnetic field vector averaged over 1/300 rotation to 
indicate the formation of Y-point around the intersection of the light cylinder and the equatorial 
plane. In a short time scales, however, we see that some part of the closed field line at the 
top of the dead zone is broken to make a magnetic island with magnetic reconnection, and 
the island go outward. However, sometimes such a island can come back and merge with the 
dead zone. Basically, plasma flows out across the outer (right) boundary. However, the plasma 
near the equator move out and in quasi-periodically. In the run, magnetic reconnection is 
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Table 2. The normalization constants and renormalization constants. 





Normalization constants 


Renormalization constants 


distance 


Xq = i?LC 


xoo = As 


time 


to = Rlc/c 


too = At/2 


velocity 


vo^c 


voo = Ax/{At/2) 


electromagnetic field 




Boo = Aa;-i/V(At/2) 


charge density 


Po = BoR^Q 


Poo = Ax-3/V(At/2) 


current density 


Jo = BoRi^qC 


Joo = Aa;-i/2(At/2)-2 


charge 


qa = BoRl^Q 


goo = Ai3/2/(Ai/2) 


mass 


mo = BlRl^c--^ 


moo = 1 



essential, and the associated heating and acceleration should be important. Figure 13 shows 
the ratio, |i^|/|-B|, averaged over the same period. The region where \E\ > \B\ degenerates 
near the equatorial plane. It is suggested that the electric-field-dominant region becomes several 
Debye length in thickness. With the initial run, we suggest that the Y-point can be active in 
magnetic field dissipation associated by magnetic reconnection. The electric-field-dominant 
region (|£^| > |-B|), which is suggested by the force- free solution, shrinks into a thin equatorial 
sheet. 

In a subsequent paper, we study the acceleration and heating, namely the conversion of 
magnetic energy into plasma quantitatively. 

d 



in 
o 
d 



O 

in 



0.995 1 1.005 1.01 

R / Rlc 

Fig. 12. The vector of the averaged poloidal magnetic field arotmd the Y-point. 
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Fig. 13. The distribution around tlic Y-point. 

4. Summary 

We develop the axisymmetric PIC code in order to analyze the Y-point in the pulsar 
magnetosphere. The code is verified by the dispersion relations of waves in electron-positron 
plasmas. We also established the appropriate boundary conditions for the Y-point in use of 
the force-free semi-analytic solution. We have performed a test run and find that the Y-point 
is actually formed at the top of the dead zone on the light cylinder. It is, however, variable 
with magnetic reconnection. We suggest that the Y-point can be a region around which the 
magnetic field energy convert into the plasma heat and outfiow and subsequently into radiation. 
We also suggest that some plasma jet can be fiow back to polar cap region to produce radio 
emission. In the next step, we will perform large scale simulations and analyze the Y-point 
structure in detail. 

Numerical computations were in part carried out on NEC SX-9 at Center for 
Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan and in 
part carried out on PRIMEPOWER 850 at Networking and Computing Service Center of 
Yamagata University. This work was supported in part by a Grant-in- Aid from the Ministry 
of Education, Science, Sports and Culture of Japan (19540235). The page charge of this paper 
is supported by CfCA. 
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